Tree method for quantum vortex dynamics 
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We present a numerical method to compute the evolution of vortex filaments in superfluid helium. 
The method is based on a tree algorithm which considerably speeds up the calculation of Biot-Savart 
integrals. We show that the computational cost scales as Nlog(N) rather than N 2 , where N is 
the number of discretization points. We test the method and its properties for a variety of vortex 
configurations, ranging from simple vortex rings to a counterflow vortex tangle, and compare results 
against the Local Induction Approximation and the exact Biot-Savart law. 
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C/5 ! I. INTRODUCTION 

o: 

Quantum turbulence [l|-|3| consists of a disordered tangle of reconnecting quantized vortex filaments, each carrying 
one quantum of circulation n. This form of turbulence is easily created by agitating superfluid liquid helium ( 4 He) 
r-j ■ with propellers 0,0], forks [f|, or grids Q, the same techniques which are used to create turbulence in ordinary fluids. 
There are also methods which are unique to superfluid liquid helium^for example heat flows Q and streams of ions 
[Tol | . Quantum turbulence is also studied in superfluid 3 He-B [TJ 121 and, more recently, in atomic Bose-Einstein 
condensate sll3l . In the study of quantum turbulence, particular attention is devoted to the similarities with ordinary 
turbulence [14l . 1 151. R ecent progress has been boosted by the development of new flow visualization techniques, such 
as tracer particles fTfjl li~7| . Andreev scattering [l8| and laser-induced fluorescence fl9j. 
' This article is concerned with numerical simulations of quantum turbulence. Star ting fro m the pioneering work of 
. Schwarz [20I ] , numerical simulations have always played an important role in the field [2lM3l| . The recent experimental 
progress, rather than reducing the need of numerical simulations, has highlighted their importance in interpreting 
^ visualization methods [H, HU . 

In superfluid 4 He, the vortex core radius (do ~ 10~ 8 cm) is many orders of magnitude smaller than the average 
separation between vortex lines (typically from 10~ 2 to 10 _4 cm) or any other relevant length scale in the flow. This 
feature was early recognised by Schwarz |20|. who proposed to model vortex lines as spaces curves s = s(£, t) of 
infinitesimal thickness, where t is the time and £ is arc length. This vortex filament approach is also valid in 3 He-B, 
although to a lesser extent because the core size is larger (ao ~ 10~ 6 cm). The space curves are numerically discretized 
by a large, variable number of points s, (i = 1, • • ■ N), which hereafter we refer to as vortex points. 

At temperatures low enough that the normal fluid and mutual friction is negligible, the governing equation of 



motion of the superfluid vortex lines is the Biot-Savart (BS) law [3 
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where the line integral extends over the entire vortex configuration C. Eq. ([1} formulates the classical Euler dynamics 
of vortex filaments in incompressible inviscid flows. The assumption of incompressibility is justified, as all relevant 
speeds, of the order of few cm/s, are much smaller than the speed of sound, which is approximately 240 m/s (but 
see the discussion of reconnections and dissipation in Section HIT)) . Unlike vortex lines in classical inviscid Euler flows, 
however, superfluid vortex lines can reconnect with each other when they come sufficiently close 35-39]. This effect is 
captured by a more microscopic model, the Gross-Pitaevskii Equation (GPE) for a Bose-Einstein condensate [40l| . In 
the context of the vortex filament approach, reconnections are modelled by supplementing Eq. (JTJ) with an algorithmical 
reconnection procedure, as originally proposed by Schwarz 20] . 

Two difficulties arise in seeking a numerical solution of Eq. ([T]). Firstly Eq. (Q]) is singular at s = r. Following 

Schwarz [2(| , to circumvent this problem we split the velocity of the filament of the i^ 1 vortex point Sj into a local 
contribution (near the singularity) and non-local contributions, writing 
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ds, k I ^JWi+\ \ , „ K I (si - r) 

-77- = -. — In s x s ' - — <h j pr x dr. (2) 

dt Air y ai J 1 1 At: J c , |s 4 - r| 3 

Here a prime denotes a derivative with respect to arc length (hence s' = ds/d^ and s" = d 2 s/d^ 2 ), and £i+i are the 
arclcngths of the curves between points S;_i and and between and Sj.fi, £ is the original vortex configuration 
C without the section between Sj_i and Sj+i, and ai is a cutoff parameter. For simplicity, hereafter we assume that 
the cutoff parameter is equal to the vortex core radius, ai ~ ao- 

To numerically calculate the nonlocal term in Eq. ([5]), we must perform further manipulations. We follow (20j 
and split the integral into discrete parts, which account for the contribution of the vortex between points along C . 
Assuming a straight line between vortex points, the contribution to the nonlocal term due to the vortex segment 
between points s,,- and s J+1 is given by, 

Aw "' (8l) = 2n(AAC-B 2 ) 

Where p = Sj — Sj, q = s J+ i — Sj, A = |p| 2 , B = 2p • q and C = |q| 2 . The total nonlocal contribution is then found 
by summing over all j € C 

The second difficulty, which we address in this paper, is the computational cost of the BS law. It is apparent from 
Eq. (JlJ that the velocity at one vortex point depends on an integral over all other N — 1 points, so the CPU time 
required to evolve the vortex configuration in time scales as N(N — 1) ~ N 2 . Because of this cost, it has been difficult 
to numerically calculate vortex tangles with vortex line density L (vortex length per unit volume) comparable to what 
is achieved in experiments. Similarly, it has been difficult to perform sufficiently long numerical calculations which 
are needed, for example, to saturate the initial vortex configuration to a steady state tangle, or to produce long time 
series suitable for spectral analysis, or to study the decay of quantum turbulence. 

Schwarz [20] proposed the Local Induction Approximation (LIA) [U, [42] as a practical alternative to the exact (but 
CPU- intensive) BS law. Neglecting the second non-local term in Eq. @, we have simply 
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where c is a constant of order unity and (R) is the mean radius of curvature of the vortex filament. The LIA greatly 
increases the efficiency of the numerical code because its computational cost scales with N, as opposed to iV 2 for 
the BS law. Moreover the LIA offers insight into vortex dynamics. Using LIA, it is easy to show that, in the first 
approximation, the motion of a point on the vortex filament is along the direction which is binormal to the filament at 
that point, with speed which is inversely proportional to the local radius of curvature R = l/|s"|. Moreover the LIA 
describes well the translational motion of a circular vortex ring and the rotation of a small amplitude Kelvin wave. 

Unfortunately, the LIA does not capture the interaction between two vortex filaments, and the interaction which 
different parts of the same vortex have on each other. For example, two circular vortex rings, set coaxially close to 
each other, move in a leapfr ogg ing fashion, but the effect is ignored by the LIA. Other limitations of the LIA have been 
discussed in the literature [43144 5j . Clearly the LIA is not the best tool to study important properties of turbulence 
which arise from the vortex interaction, such as how the energy is distributed over the length scales (the Kolmogorov 
spectrum). 

In a recent paper [3l| we presented the spectrum of vortex density fluctuations in superfluid turbulence and compared 
it to experiments. This calculation was possible because we used a tree algorithm, which retains the long-range 
interaction of the BS law, but sensibly approximate it, so that its computational cost scaled as TV log (N) rather than 
N 2 . The difference between N 2 and Nlog(N) cannot be overstated: it is the same improvement of the Fast Fourier 
Transform. For example, in the work reported in ref. [3l| . the typical number of vortex points which we used was of 
the order of N ~ 10 5 . Tree algorithms are popular in astrophysical N-body simulations [461 ]. The numerical solution 
of the gravitational interaction of N bodies, in fact, has the same difficulty of the BS law: the motion of one body 
depends on the other N — 1 bodies, hence the computational cost scales as N 2 , as for vortex points. 

The aim of this paper is to present a detailed description of a tree algorithm tailored to vortex dynamics, provide 
quantitative error estimates associated with the method, and compare results with methods used in the literature. 
In particular we show that the tree method provides huge computational savings but without the inherent errors 
associated with the LIA. 

The structure of the article is as follows. Section [II] outlines the algorithms and methods required to create a tree 
structure for a set of points embedded in three-dimensional space. We pay careful attention to how this tree structure 
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FIG. 1: (Color online) Illustration of the tree construction in two dimensions (quad-tree). The points (blue dots) are enclosed 
in the root cell (a), which is divided into four cells of half size (b), until (c,d) there is only one point per cell. 



can be efficiently used to compute the velocity at a particular vortex point. In Section IlIII we provide further details 
on the algorithms required to apply the vortex filament method to the study of quantised turbulence. Section IIVI 
contains the results of numerical tests of the tree algorithm. We show that the N 2 scaling inherent in the evaluation of 
the BS integral can be approximated, with a tolerable loss of accuracy, by a tree method which exhibits an improved 
NlogN scaling. A recent, and very detailed, numerical study of counterflow turbulence by Adachi et al. [45[ as again 
demonstrated that the LIA is unsuitable for the study of quantum turbulence. In Section [V] we verify that the tree 
method does not exhibit the same deficiencies of the LIA. 

In all simulations presented, here we use parameters which refer to superfluid 4 He: circulation k = 9.97 x cm 2 /s 
and vortex core radius ao ~ 10~ 8 cm. All our results and methods, however, can be generalised to vortex systems in 
low temperature 3 He-B. 



II. TREE ALGORITHM 

Building on the pioneering work of Barnes and Hut (49| , the majority of astrophysical and cosmological N-body 
simulations have made use of tree algorithms to enhance the efficiency of the simulation with a relatively small 
loss in accuracy [Hfj. The advantages of these methods are two- fold. Firstly they have been shown to scale as 
iVlog(iV)), where N is the number of particles used in the simulation. Secondly there are now openly available 
parallel implementations of the tree algorithm, which have been used to compute simulations with up to 10 billion 
particles [5l|. Our implementation of the tree algorithm for vortex dynamics is limited to a serial program, however 
the development of a parallel scheme is one of our future aims. 



We are not the first to apply tree methods to quantised vortex simulations. Kivotedes [47] used a tree method to 
enhance a vortex filament simulation, but did not provide any description, testing and evaluation of the performance 
of the method. Kozik and Svistunov 48] introduced an interesting method to speed up the simulation of the Kelvin 
wave cascade along quantised vortices using a scale separation scheme. Whilst this is a significant improvement over 
the full BS integral the problem discussed in |48| is inherently one-dimensional and it is not clear if there is an efficient 
implementation of the scheme that would enhance the speed of a three dimensional simulation. 

Detailed descriptions of the tree algorithm can be found in a number of works (52l . [54|, but their focus is on 
astrophysical simulations. In this paper we provide a detailed description of the method which is accessible to the 
fluid dynamics and the quantum turbulence community. 

At each time-step we have a set of orientated vortex points which represent the current position and direction of the 
vortex filaments. We first group the vortex points into a hierarchy of cubes which is arranged in a three-dimensional 
octree structure. The octree can be either be constructed top down or bottom up. We think that it is conceptually 
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FIG. 2: The complex octree structure required to map the vortex points for a single vortex loop (left panel) is plotted in the 
middle panel. The right panel shows the same structure in the :ri/-plane,at z = 0. 




FIG. 3: A cartoon of the calculation of the centre of circulation (green) triangle and the total circulation bold gray arrow. The 
(red) circles mark the position of the vortex points within a tree cell and the black line gives a representation of the vortex 
filament. 



easier to work from the top down, repeatedly dividing the domain. To achieve this aim, we first divide the full 
computational box, the root cell, into eight cubes, and then continue to divide each cube into eight smaller 'children' 
cubes, until either a cube is empty or contains only one vortex point. 

The numerical code, QvORT[63j we have developed is written in FORTRAN 2003 and employs linked lists [53| 
to efficiently create the tree structure. In the code each cell in the tree is a node which, using pointers, links to 
the 8 children within the domain of the parent node. The main node is set to the full computational box, we then 
repeatedly call a recursive subroutine to allocate the points in this box to the main nodes child cells. This process is 
then repeated dividing the main nodes child cells until the hierarchical tree structure is created. 

As we create the tree, we calculate the total circulation contained within each cube and the 'center of circulation' 
using the points that the cube contains. We denote the number of points the cube contains as N c , then the total 
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circulation is given by 

No 

s s =J2(s j+1 -Sj). (5) 

3=1 

Likewise the center of circulation is calculated as, 

1 Nc 

3=1 

A schematic diagram for these quantities is displayed in Fig [31 

The time required for constructing the tree scales as Nlog(N), so it feasible (and essential) to 'redraw' the tree at 
each time-step. Fig. [T] illustrates this procedure in two dimensions for clarity. Fig. [5] shows the complex nature of the 
octree required for a single vortex ring in a three-dimensional calculation. 

The enhanced speed of the tree method appears when we calculate the induced velocity at each vortex point Si. 
As opposed to the BS method, we only take the full contribution of points near to the vortex point in question. The 
induced velocity from far vortex points is an average contribution, hence the number of evaluations required per point 
is significantly smaller than N — 1. Therefore for each point we must 'walk' the tree, and decide if a cell in the octree 
is sufficiently far. 

This is done using the concept of the opening angle 9. In their original work [49], Barnes & Hut calculated 9 by 
evaluating the distance of the center of vorticity from Sj, which we denote as d. Since each cube in the octree has 
a known width w, then the opening angle is 9 = w/d. We then must define a critical (or maximum) opening angle 
9max- If 9 < 6* ma x we accept the cube, and its contribution is used in computing the nonlocal contribution to the 
velocity via Eq. ([3]). The crucial difference between the full BS law and the tree approximation is that the vectors p 
and q must be redefined as p = s — s, and q = ss- If 9 max — then each cell contains one point at the most, and it 
is clear that = s J+ i — Sj, s = sj and we recover p and q used in the full BS integral as one would expect. 

If 9 > 9 maxi then we open the cube (assuming it contains more that one vortex point) and repeat the test on each of 
the child cubes that it contains. As a link list is used, from a programming point of view, it is relatively easy to start 
at the main node and follow the separate branches of the tree, at each stage evaluating whether the cells contribution 
can be used. The tree-walk ends when the contributions of all cubes have been evaluated. We note that the adjacent 
points are exempt from this walk, and their contribution is taken as in Eq. ^ so as to avoid the singularity. 

A correction to this method was proposed by Barnes (5^. [55j to avoid errors which would arise if the center of 
vorticity is near the edge of the cube. This requires the calculation of the distance from the center of vorticity and 
the geometrical center of the cube, which we denote as (. The test criteria then becomes 9 = w/(d — (). In Section 
IIVI we evaluate the improvement caused by this new criterion. 

Whatever opening criterion is used, it is crucial to have a reasonable choice for #max- If $max is too large, the 
discrepancy between the full BS integral and the tree approximation will damage to the accuracy of the simulation. 
In Section IIVI we shall also show how the errors in the calculation of the velocity at the vortex points depends on 
9m&x- 

In the next section we describe other necessary details of the vortex filament method which we used. 



III. NUMERICAL METHODS 

The numerical simulations presented here use a 3 rd order Adams-Bashforth scheme. Given an evolution equation 
of the form ds/dt — v, the recursion formula is 

si 14-1 = $ + ^(23< - lGvr 1 + 5v»- 2 ) + O(At^). (7) 

where At is the time-step and the superscript n refers to the time t n = nAt (n = 0, 1, 2, • • •). Lower order schemes 
are used for the initial steps of the calculations, when older velocity values are not available. 

As the vortex points move, the separation between neighbours along the same filament is not constant. In general 
this is not a problem, as we can use finite-difference schemes which work with an adaptive mesh size. However 
the distance between points sets the resolution of the calculation, therefore we must set some upper-bound on the 
distance between vortex points before we re-mesh the filaments by introducing new vortex points. Our criterion is the 
following: if the separation between two vortex points, and s^+i, becomes greater than some threshold quantity S 
(which we call the minimum resolution), we introduce a new vortex point at position Sii given by 
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Si< = \(si + s i+1 ) + Ub*. Ri,j , (8) 

where B41 — js^l -1 . In this way s ir = (s i + s i+1 )/2, that is to say the insertion of new vortex points preserves the 
curvature. This property is desirable, as the curvature of the filament directly affects the velocity field via Eq. ((2J. 
Our procedure is thus superior to simpler linear interpolation. 

Conversely, vortex points are removed if their separation becomes less than 6/2, ensuring that the shortest length- 
scale of the calculation remains fixed. This property is important, as the maximum time-step which can be used is 
related to this length-scale. In the long-wavelength approximation (fcao <C 1), the angular frequency of a Kelvin wave, 
along a straight vortex is 

W «~Pn2/(fcoo)- 7 ], (9) 

where 7 = 0.5772 is Eulers constant. Hence the timescale for the fastest motions, with a maximum wavenumber 
k max = Att/6, is 

x 1 A5/2) 2 , . 

W 1 ~ , 1 (hi v 10 

4Kln(d/27rao — 7) 
thus a reasonable maximum time-step that can be used is given by 

M< P/f v (11) 
Klog(J/27ra ) V ; 

We note that at finite temperatures this condition may be relaxed as mutual friction will tend to dampen Kelvin waves. 
Setting a consistent minimum separation between vortex points makes the numerical resolution of the calculation 
transparent, and prevents the creation of too short length-scales which would considerably slow down the calculation. 

In order to model any turbulent system we must pay particular attention to the dissipation of kinetic energy. Even 
at very low temperatures, in the absence of mutual friction with the normal fluid, there are still mechanisms which 
dissipate the kinetic energy of quantised vortices: reconnections, decay of very small vortex loops and direct phonon 
emission. 



Using the GPE condensate model, Leadbeater & al |58[ showed that, at vortex reconnection events, a small amount 
of kinetic energy is transformed into acoustic energy of a rarefaction pulse which moves away. Since our model neglects 
sound waves and it would be impractical to determine with accuracy very small changes of kinetic energy, we interpret 
this dissipation of kinetic energy as a reduction of vortex length, essentially using length as a proxy for energy. 

We model reconnection events algorithmically in the following way. If two discretization points, which are not 
neighbours, become closer to each other than the local discretization distance, a numerical algorithm reconnects 
the two filaments by simply switching flags for the vortex points in front and behind the filament, subject to the 
criterion that the total length decreases. Self-reconnections (which can arise if a vortex filament has twisted and 
coiled by a large amount) are treated in the same way. Since reconnections involve only anti-parallel filaments, prior 
to reconnection we form local (unit) tangent vectors s', and, using the inner product, we check that the two filaments 
are not parallel. 

Clearly, testing for possible reconnections requires a distance evaluation between a vortex point and the vortex 
points in the vicinity. The most efficient way to do this is within the tree- walk itself, creating a linked list of all vortex 
points within a radius of 6/2 from a vortex point in question. This avoids the need to re-loop over all discretization 
points, which would involve approximately N 2 operations. 

We distinguish two kinds of phonon emissions. The condensate model shows that very small vortex loops, of radius 
of the order of do, which move at speed close to the speed of sound, are unstable and decay into sound 59]. It would 
not be practical to mesh vortex filaments down to this atomic scale, so we remove any small vortex loop with less 
than five discretization points. 

Phonon emission can also be induced by very short Kelvin waves which rotate rapidly [60L [6l| . Again, since we 
do not mesh filaments down to this scale, we remove vortex points if the local wavelength is smaller than a specified 
value [3l[. Since our smallest length-scale is maintained at 6/2, the maximum curvature is of the order of 2/6. We 
proceed in this way: at each time-step we compute the local curvature C(£) = |s"| at each vortex point. If, at some 
location along a vortex filament, the local curvature exceeds the critical value 1.9/6 cm -1 (95% of the maximum value 
2/5), this vortex point is removed, smoothing the vortex filament locally. This loss of length as a proxy of energy is 
our model of phonon emission in the low temperature limit. 
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We approximate all spatial derivatives using 4 th order finite difference schemes which account for varying mesh sizes 
along the vortex filaments 62]. A high-order scheme is desirable as the spatial derivatives enter into the equation 
of motion via, Eq. ([2]). Let Sj be the i th point on the vortex filament; the two points behind have positions s,;_2 
and Sj_i, and the two points in front have positions s$+i and Sj+2. We denote 4-1 = |sj_i — Sj—al, 4 = |sj — s»_i|, 
4+1 = \si+i - Si|, and 4+ 2 = \s i+2 - s l+1 \. We have: 

s- = AuSi^ 2 + BuSi-i + CijSj + Dus i+ i + E u s l+2 , (12) 

where the coefficients An, Bn, Cu, Du and En are given by, 

. = 44^+44+14+2 

11 4-l(4-l+4)(4-l+4 + 4+l)(4-l+4 + 4+l+4+ 2 ) 1 ' 

„ _ — 4-l4+l ~ 44+1 ~ 4— l4+l4+2 — 44+l4+2 , 
U ~ 4-l4(4+4+l)(4 + 4+l+4+ 2 ) ( ' 

Cxi = -(Au + Bxi + D u + Exi) (15) 

D _ 4-l44+l + g4+l + 4-l44+2 + gjlj+2 , . 

" 4+l4+ 2 (4 + 4+l)(4-l+4+4+l) [ ' 

E = -h+lij - 4-l44+l , , 

11 4+2(4+l+4+2)(4+4+l+4+2)(4-l+4+4+l+4+ 2 ) 1 ' 

Similarly we have: 

s" = A 2i Si^2 + B 2 iSi-i + C 2l s l + D 2 iS l+ i + E 2 iS t+2 , (18) 
where the coefficients An, B2i, C 2 i, D 2 i and E2i are given by 

= 2[-2iji i+ i + - 44+2 + Mi+j 

21 ~4-i(4-i+4)(4-i + 4 + 4+i)(4-i +4 + 4+1 +4+2) ( ' 
_ 2[24-i4+i + 244+i - + 4-i4+2 + 44+2 - W . , 

21 _ 4-i4(4 + 4+i)(4 + 4+i+4+ 2 ) 1 j 

C 2i = -(^2i + B 2i + .D 2i + E 2i ) (21) 

= 2[-4-i4 - if + 4-i4+i + 24-4+1 + 4-14+2 + 24-4+2] f , 

22 4+i4+ 2 (4 + 4+i)(4-i +4 + 4+i) 1 j 

B = 2[4-i4 + ^-4-i4+i -244+1] , . 

4+2(4+1 +4+2) (4 + £ i+i + 4+ 2 )(4-i + 4 + 4+i + 4+2) 

Note that if we set 4-1 = 4 = 4+1 = 4+2 = then the above expressions reduce to familiar finite-difference 
schemes for a uniform mesh: 

= y^( s »-2 - 8s,_i + 8si+i - s l+2 ) + Oih*), (24) 

s'/ = + 16si_i - 30s 4 + 16s i+1 - s l+2 ) + 0{h A ). (25) 

The final issues to address in this section is the boundary conditions used in the simulations and how they are 
implemented. All numerical simulations presented here were performed with periodic boundary conditions. The 
computational domain was a cube. If a vortex point leaves the domain, then it is simply re-inserted at the opposite 
side of the cube. The presence of periodic boundary conditions affects the calculation of the distance between two 
vortex points (two vortex points point at the opposite end of the domain may be actually close to each other). We 
must also use periodic 'wrapping' to ensure that the velocity field is periodic; in theory we should wrap the cube 
with an infinite number of copies in all directions. In practice, we only apply a single layer of wrapping, which still 
requires a very large amount of numerical evaluations. If we used the BS law, a single layer of wrapping would 
require 27N 2 evaluations per time-step. Fortunately again the tree method eases the amount of numerical work, as 
the computational domain is wrapped with copies of the octree we have already computed. Walking each of these 
trees, as well as the tree for the main domain, does not require a significant increase of CPU time. 
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FIG. 4: The CPU time required for one hundred time-steps, plotted as a function of the number of vortex points N. The 
calculations were performed on an single processor with a clock speed of 2.53 GHz. The symbols represent different values of 
flmax as follows: triangles, #max = 0.6, asterisks, #max = 0.4, squares, #max = 0.2, and circles, #max = 0. (that is, the full 
BS calculation). 



IV. TESTS OF THE TREE ALGORITHM 



We now test the error associated with the use of the tree approximation, and show that, with an appropriate choice 
of #max, the discrepancy between using the tree algorithm and the full BS law is very small. All tests are carried out 
in a periodic cube of side 0.1 cm. The minimum resolution is set at 6 = 0.001cm, and the time-step is At = 10 _5 s. 
For our first test we consider a system of ten random vortex loops with an increasing number of vortex points (hence 
increasing radius). We determine the CPU time required for evolving the system for one hundred time-steps as a 
function of 9; the results are presented in Fig. [4] It is apparent that even for very small maximum opening angles, 
#max, there is a huge increase in CPU speed, and for larger values of 8 m a,x the A log N scaling is apparent. 

The second test consists in determining the discrepancies between the 'true' velocity field (obtained using the BS 
law) and the velocity filed obtained with the tree method, as a function of 6* m ax- Let v^g be the the velocity at 
the discretization point Sj obtained using the BS integral, and v^x the corresponding velocity obtained using the tree 
method. We define the mean relative percentage error of the tree method as 



100 x 



|v ijBS | 



(26) 



Again, we use a system of 10 randomly oriented and positioned vortex loops in the same periodic box, each vortex 
loop discretized into 200 vortex points. Figure [5] shows how e varies as a function of #max, for both the original tree 
opening criterion and the correction of Barnes. Whilst the improvement is modest, since the extra computational 
cost is negligible, hereafter we include Barnes's modification. 

The next test investigates how this error varies as a function of N for a fixed value of max - This test is particularly 
important in view of computing very dense vortex tangles, which are not accessible using the exact BS law. As before, 
we consider a system of 10 random vortex loops with an increasing number of vortex points in the same periodic box, 
and compute e. Figure [6] shows the results for fixed #max = 0.4. It is apparent that the error grows with the number 
of vortex points, but its magnitude is still small considering the advantages of the tree method. We obtain the same 
result with a different set up to increase N: we fix the vortex loop size but increase the number of vortex loops. 

In the previous tests we have only considered the errors at a single time-step of a simulation. It is important to 
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#max 



FIG. 5: The scaling of the mean relative percentage error in the velocity field, e, evaluated at discretization points plotted 
as a function of #max- The dashed line corresponds to the original opening angle criterion [49[ and the solid line shows the 
performance of the correction of Barnes [H^] . 
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FIG. 6: The mean relative percentage error of the tree approximation for a fixed maximum opening angle, #max = 0.4, with 
an increasing number of discretization points and therefore vortex length in the system. 
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FIG. 7: A system of leap-frogging vortices, each loop contains 200 vortex points with S — 0.001cm, the mean separation is 
0.00075cm giving each loop a radius of 0.024 cm. These snapshots are taken from a tree code simulation with #max = 0.4. 



quantify how the difference between the exact BS law and the tree method varies in time. For this test we deliberately 
choose a vortex configuration which the LIA fails to accurately model, namely two 'leap-frogging' vortices. 

The initial condition consists of two co-axial vortex rings set in the xy-plane, separated by a small distance in the z 
direction. Given the orientation of the vortex lines, the local velocity component - see Eq. ^ - drives the vortices in 
the negative z direction. The non-local velocity component makes the vortex ring which is initially behind to squeeze 
through the vortex ring which is in front, overtaking it; the process repeats over and over again, in a leap-frogging 
fashion (note that the vortices maintain the circular shape) . This remarkable behaviour is visible in Fig. [7J We stress 
that this behaviour is not captured by the LIA. We repeat the calculation using the exact BS law, keeping the same 
number N of vortex points and the same minimum resolution S. We then quantity the difference between the tree 
method and the BS law by computing 



D,W = ±± KB3{t) ; Bi ' m , (27) 

2 — 1 

where s^ss is the position of the fi 1 vortex point in the BS calculation, and s^t is the corresponding position in 
the tree method calculation. The temporal duration of these calculations captures four separate leap-frogging events. 
Figure [8] shows the evolution of D s for two different values of #nmx: 0.4 and 0.8. It is apparent that in both cases 
the difference between the tree method and the exact BS evolution is small (smaller than the minimum discretization 
level <$), particularly for #max = 0.4. 



V. COUNTERFLOW SIMULATIONS 



Having benchmarked the tree method for simple vortex configurations, we change our focus to numerical simulations 
of quantum turbulence. We concentrate our attention to turbulence sustained by the relative motion (countcrflow) of 
normal fluid and superfluid components, driven by a heat flow at non nonzero temperature. This form of turbulence 
which has no analogy in ordinary fluids was originally studied experimentally by Vinen 8], and, numerically, by 
Schwarz [201 ] . Counterflow turbulence has been recently re-examined by Adachi & al. [451 ] . who compared results 
obtained using the LIA (as in the pioneering work of Schwarz 20]) and the full BS law. Counterflow turbulence is 
particularly suitable for testing the tree method, as Adachi & al. used it to discuss the shortcoming of LIA. We want 
to make sure that the tree method does not suffer from the same problems of the LIA. 

At finite temperatures we must account for the mutual friction between the superfluid and the normal fluid. The 
equation of motion at position s is given by 

ds 

— = v s + as' x (v„ - v a ) - a's' x ]s' x (v n - v s )] , (28) 
at 

where a, a' are temperature dependent friction coefficients [56l[57j. v„ is the normal fluid's velocity, and the velocity 
v s is calculated using the tree approximation to the de-singularised BS integral. As usual, the evolution includes the 
algorithmical vortex reconnection procedure. 
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FIG. 8: The evolution of the normalised distance error statistic D s (see Eq. ((27|), for a system of leap-frogging vortex rings, 
plotted as a function of t (seconds). The solid line corresponds to a simulation with (9max = 0.4, and the dashed line to 
#max = 0.8. 



The calculation is performed in a periodic cube with sides of length 0.1 cm as done by Adachi et al. [45]. Superfluid 
and normal fluid velocities v„ and v s are imposed in the positive and negative x directions respectively, where 
v n s = | v n — v s | is proportional to the applied heat flux. Simulations are performed with the same numerical resolution 
6 = 8 x 10~ 4 cm and time-step At = 10~ 4 s as used by Adachi et al. at various value of v n and T. Whereas Adachi et 
al. used the BS law, we use our tree method with #max = 0.4. The number of vortex points used in these calculations 
is not large, in order to overlap with calculations which used the BS law:) for example a typical value is N rj 6000 at 
v ns = 0.572 cm/s and T = 1.9 K. It must be stressed that, as in the works of Schwarz and Adachi et al. , the back- 
reaction of the superfluid vortex lines onto the normal fluid velocity is neglected (an effect which may be significant 
at low temperatures). The initial condition consists of six vortex loops set at random locations, as for Adachi et al. . 
We find that the results do not depend on the initial condition. 

During the evolution we monitor the vortex line density 

*> = ?f c #> ( 29 ) 

where V is the volume of the computational domain and C is the entire vortex configuration. 

We find that, after an initial transient, the vortex line density saturates to a statistical steady state value, see 
Figure IH1 (left), which agrees well with the value reported by Adachi et al. in their Figure 2 at the same temperature 
and counterflow velocity. Figure (right) shows the expected relation \f~L ~ v ns between the steady state value of L 
and the driving counterflow. More precisely, we obtain 7 = 140.8, which compares well with 7 = 140.1 obtained by 
Adachi et al. , where L = 7 2 v^ s . 

Another important quantity considered by Schwarz and Adachi et al. is the anisotropy parameter 



Z„ = ±.^[l-( s '.r {{ ) 2 ]d4, (30) 

where fn is the unit vector parallel in the direction of the normal fluid. Figure [T0l (right) shows that Z11 is almost 
independent of v ns , in agreement with experiments and with the numerical calculations of Adachi et al. performed 
using the full BS law (shown in their Fig. 7). For example at v ns = 0.55 cm/s we have I\\ ~ 0.83, compared to 0.80 - 
0.85 of Adachi et al. . Figure ITD1 also shows that this anisotropy is not captured by the LIA. 
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FIG. 9: Left: The evolution of the vortex line density L (cm -2 ) vs time t (s) for counterflow simulations at 1.9 K for four 
counterflow velocities: dashed-dotted, v na = 0.2 cm/s, dashed, v n s = 0.4 cm/s, solid (black), v ns — 0.6 cm/s, and solid (grey), 
Vns = 0.8 cm/s. Right: The plot of the vortex line density L (cm -2 ) vs v ns (cm/s) confirms that the saturated line density 
scales according to L ~ v 2 s . 
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FIG. 10: Left: The anisotropy parameter lu plotted as a function of time t (s) for the LIA simulation (grey) and for the tree 
method (black) with 6*max = 0.4. Right: The steady-state hi parameter plotted for four different normal fluid velocities v na 
(cm/s) plotted with error bars (calculated from temporal fluctuations) for the tree method only. These results are in excellent 
agreement with those of Adachi et al. Hfl 



A final qualitative test is shown in Fig. 1111 the left panel shows the bundling of vortex lines wrongly predicted by 
the LIA, and the right panel shows the vortex tangle computed using the tree method, in agreement with Adachi et 
al. . 



VI. CONCLUSIONS 



We conclude that the tree method is an ideal tool to perform numerical simulations of quantum turbulence. It does 
not suffer from the known shortcoming of the LIA, retaining to a good approximation the nonlocal interaction which 
is typical of the exact BS law. Unlike the BS law, whose computational cost scales as A^ 2 where N is the number 
of discretization points, the cost of tree method scales as N\og(N), thus making possible numerical calculations of 
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FIG. 11: Views of the tangles from the top of the computational domain in the xy-plane for the LIA simulation (left) and the 
tree method (right). The domain is a cube with sides of length 0.2 cm; the imposed normal velocity is v na = 0.5 cm/s and 
T = 1.6 K. Note the stratified nature of the tangle wrongly predicted by the LIA. 



intense vortex tangles for a sufficiently long time. Finally, the tree method can be further speeded up by parallelization. 
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